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ABSTRACT 


Accurate  seismic  event  location  is  integral  to  the  effective  monitoring  of  the  Comprehensive  Nuclear-Test-Ban 
Treaty  (CTBT),  as  well  as  being  a  fundamental  component  of  earthquake  source  characterization.  To  account  for  the 
effects  of  crustal  and  mantle  structure  on  seismic  travel  times,  and  to  improve  seismic  event  location  in  the  Middle 
East  and  North  Africa  (MENA),  we  are  developing  a  set  of  radially  heterogeneous  and  azimuthally  invariant  travel¬ 
time  models  of  the  crust  and  upper  mantle  for  each  MENA  seismic  station. 

We  begin  by  developing  an  average  one-dimensional  velocity  model  that  minimizes  the  P-phase  travel-time 
residuals  from  regional  through  teleseismic  distance  at  each  station.  To  do  this  we  (1)  generate  a  suite  of  1-D 
velocity  models  of  the  earth,  (2)  compute  travel  times  through  the  1-D  models  using  a  tau-p  formulation  to  produce 
standard  travel-time  tables,  and  (3)  minimize  the  root-mean-square  (rms)  residuals  between  the  P-phase  arrivals 
predicted  by  each  model  and  a  groomed  set  of  ISC  P-phase  arrival  times  (Engdahl  et  al.,  1998).  Once  we  have  an 
average  one-dimensional  velocity  model  that  minimizes  the  P-phase  travel-time  residuals  for  all  distances,  we  repeat 
steps  1  through  3,  systematically  perturbing  the  travel-time  model  and  using  a  grid  search  procedure  to  optimize 
models  within  regional,  upper  mantle,  and  teleseismic  distance  ranges.  Regionalized  models  are  combined  into  one 
two-dimensional  model,  using  indicator  functions  and  smoother  methodologies  to  reduce  distance  and  depth 
discontinuity  artifacts  between  the  individual  models. 

Preliminary  results  of  this  study  at  a  subset  of  MENA  stations  show  that  we  are  improving  predictability  with  these 
models.  Cross-validating  the  travel-time  predictions  with  an  independent  data  set  demonstrates  a  marked  reduction 
in  the  variance  of  the  travel-time  model  error  distributions.  We  demonstrate  the  improvement  provided  by  these  2-D 
models  by  relocating  the  1991  Racha  aftershock  sequence.  We  will  extend  our  investigation  to  additional  MENA 
stations,  and  will  use  our  model  in  tandem  with  nonstationary  empirical  corrections  (nonstationary  Bayesian  kriging) 
to  further  improve  our  ability  to  accurately  predict  travel  times  and  locate  seismic  events  in  this  region. 
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OBJECTIVE 


An  accurate  velocity  model  of  the  earth  is  a  fundamental  component  of  seismic  event  location.  Global  velocity 
models  such  as  akl35  (Kennett  el  al.,  1995)  minimize  a  set  of  worldwide  travel-time  residuals  and  are  primarily 
meant  for  prediction  of  teleseismic  travel-times.  However,  global  models  are  often  inadequate  for  prediction  of 
travel  times  at  regional  and  near  teleseismic  distances,  where  the  effect  of  crustal  and  mantle  structure  on  seismic 
travel  times  is  considerable.  In  order  to  improve  the  travel-time  predictions  at  regional  and  near  teleseismic 
distances  in  the  Middle  East  and  North  Africa  (MENA),  we  are  developing  two-dimensional,  station-specific 
velocity  models  that  are  optimized  to  predict  travel  times.  We  will  use  our  travel-time  models  at  each  station  in 
tandem  with  nonstationary  spatial  corrections  (nonstationary  Bayesian  kriging)  to  further  improve  our  capability  to 
accurately  locate  all  seismic  events  in  this  region. 
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RESEARCH  ACCOMPLISHED 


Data 

We  focus  our  initial  travel-time  modeling  efforts  on  27  of  the  International  Monitoring  System  (IMS)  seismic 
stations  in  the  Middle  East  and  North  Africa.  The  locations  of  the  IMS  stations  used  in  this  study  are  shown  in 
Figure  1 .  The  data  we  use  are  a  set  of  E-phase  arrival  times  from  earthquakes  recorded  at  those  each  of  those 
stations.  The  earthquakes  are  a  relocated  subset  of  those  contained  in  the  ISC  Bulletin  (Engdahl  el  al.,  1998).  To 
minimize  the  error  introduced  by  a  misidentification  of  phase  arrival  time,  we  use  only  high-confidence  E-phase 
picks.  The  data  set  has  also  been  declustered  -  a  statistical  grooming  procedure  that  reduces  the  impact  of  outliers, 
enhances  numerical  stability,  and  lessens  computation  demands.  We  divide  the  groomed  data  set  into  a  modeling 
data  set  and  a  cross-validation  data  set.  The  cross-validation  data  set  represents  10%  of  the  total  data  set. 


Method 

For  the  purposes  of  this  study,  the  data  are  parsed  into  three  earthquake-station  distance  ranges:  regional  (1°-13°), 
upper  mantle  (13o-30°),  and  teleseismic  (30o-90°).  Sensitivity  analyses  were  performed  to  identify  those  properties 
of  the  crust  and  mantle  to  which  seismic  travel  times  in  each  distance  range  were  most  influenced.  Of  the  1 1 
properties  investigated  in  the  sensitivity  analyses,  crustal  thickness,  upper  and  lower  crustal  E- wave  velocity  and 
upper  mantle  velocity  had  the  largest  effect  on  travel-times.  The  eight  model  parameters  we  use  to  describe  model 
space  are  listed  below  in  Table  1.  We  develop  an  adaptive  grid  search  method  that  efficiently  samples  the  space  of 
reasonable  models,  allowing  the  four  most  influential  model  parameters  (as  identified  in  the  sensitivity  analyses)  to 
vary. 


Table  1:  Grid  search  model  parameters  of  the  1-D  velocity  models 


MODEL  PARAMETERS  VARIED  IN  THE  GRID  SEARCH 
Crustal  thickness  30-55  km 
Crustal  E- wave  velocity 

Upper  crust  5. 5-7.0  km/s 
Lower  crust  6. 6-8.0  km/s 
Upper  mantle  velocity  7.9-8. 1  km/s 

ADDITIONAL  MODEL  PARAMETERS 
Sediment  thickness  4  km 
Sediment  E- wave  velocity  4  km/s 
Thickness  of  mantle  lid  25  km 
E- wave  velocity  gradient  in  the  mantle  lid  0 


We  calculate  travel  times  through  each  regionalized  E- wave  velocity  model  using  a  ray-tracer  that  employs  the 
single-valued  tau-p  formulation  similar  to  that  of  Buland  and  Chapman  (1983).  An  earth-flattening  transformation 
is  used  to  account  for  the  sphericity  of  the  earth,  which  preserves  the  kinematic  properties  of  the  rays.  The  resulting 
travel-time  tables  are  populated  with  travel  times,  parameterized  by  distance  and  depth. 

To  test  the  predictive  power  of  each  model,  we  compute  the  root-mean-square  and  mean  residuals  between  the 
declustered  E-phase  arrivals  at  each  station  and  the  arrivals  predicted  by  each  1-D  regionalized  model.  We  repeat 
the  calculation  for  the  declustered  E-phase  arrivals  at  each  station  and  the  arrivals  predicted  by  the  akl35  model 
(Kennett  el  al.,  1995).  To  compute  these  residuals  we  interpolate  between  grid  nodes  to  calculate  the  predicted 
travel  time  for  the  earthquake-station  path.  Then  each  predicted  time  is  subtracted  from  the  observed  arrival  time. 
We  find  that  a  suite  of  models  predict  the  travel-time  arrivals  equally  well.  To  further  optimize  our  models,  the  50 
models  that  provide  the  lowest  rms  residuals  are  further  minimized  against  the  akl35  model.  The  model  that  most 
closely  approximates  ak!35  is  designated  as  the  preferred  model. 


The  preferred  one-dimensional  regional  models  are  then  merged  into  one  two-dimensional  model  over  the  entire 
distance  range  of  the  model,  from  regional  out  to  teleseismic  distances  (Figure  2).  Merging  is  accomplished  using 
indicator  functions  to  reduce  distance  and  depth  discontinuity  artifacts  between  the  individual  models  (Figure  2). 
The  result  is  a  radially  heterogeneous  and  azimuthally  invariant  travel-time  model  of  both  the  crust  and  upper 
mantle.  This  methodology  provides  optimal  models  for  the  three  distinct  ray-bottoming  depths,  allowing  increased 
predictability  and  a  smooth  travel-time  curve. 

To  test  the  predictive  power  of  the  2-D  models,  we  compute  the  rms  and  mean  residuals  between  the  declustered  P- 
phase  arrivals  at  each  station  and  the  arrivals  predicted  by  both  the  optimized  2-D  model  and  the  akl35  model. 
Station-specific  analyses  of  the  capability  of  our  2-D  models  to  predict  the  observations  help  us  identify  regions 
where  the  model  needs  to  be  improved  or  alternative  models  need  to  be  applied. 

There  are  several  caveats  associated  with  our  technique  that  merit  mention.  First,  since  the  travel-time  modeling 
process  is  non-linear  and  the  number  of  model  parameters  is  limited,  the  solution  is  non-unique.  In  addition,  these 
two-dimensional  models  are  azimuthally  invariant  and  radially  segmented  to  facilitate  statistical  averaging.  As  a 
result,  the  models  do  not  account  for  azimuthal  changes  in  structure  and  can  only  account  for  average  changes  in 
radial  structure.  Our  two-dimensional  models  are  thus  likely  to  have  better  predictive  power  in  areas  with  dense  GT 
event  coverage.  Stations  with  small  azimuthal  deviations  benefit  most  from  these  corrections.  To  predict 
corrections  in  aseismic  regions,  three-dimensional  models  will  likely  be  required. 


Results 

For  the  majority  of  the  27  IMS  stations  for  which  we  have  computed  2-D  models,  our  2-D  model  predicts  the 
observations  very  well,  showing  significant  variance  reduction;  for  a  small  number  of  stations,  the  2-D  model  only 
slightly  reduces  the  variance.  We  have  selected  four  IMS  MENA  stations  to  illustrate  in  more  detail  our  method  and 
results:  AAE  (Addis  Ababa,  Egypt),  ANTO  (Ankara,  Turkey),  AQU  (Abruzzo,  Italy)  and  TBT  (Canarias,  Spain) 
(Figure  1).  The  results  from  these  individual  stations  are  fairly  representative  of  results  from  the  remainder  of 
stations  included  in  our  initial  modeling  study. 

The  map  in  the  upper-right-hand  corner  of  Figure  3  shows  the  azimuthal  distribution  of  a  typical  modeling  data  set. 
This  particular  data  set  was  recorded  at  station  AAE.  The  data  have  been  parsed  into  three  earthquake-station 
distance  ranges:  regional  (white  circles;  1°-13°),  upper  mantle  (gray  circles;  13°-30°),  and  teleseismic  (black  circles; 
30°-90°). 

We  tested  the  predictive  power  of  each  P- wave  velocity  model  at  stations  AAE ,  ANTO ,  AQU  and  TBT  by  computing 
the  rms  and  mean  residuals  between  the  declustered  P-phase  arrivals  at  each  station  and  the  arrivals  predicted  by 
each  1-D  regionalized  model.  We  repeat  the  calculation  for  the  declustered  P-phase  arrivals  at  each  station  and  the 
arrivals  predicted  by  the  akl35  model  (Kennett  et  al.,  1995).  The  50  models  that  provide  the  lowest  rms  residuals 
are  further  minimized  against  the  akl35  model.  For  each  station,  the  model  that  most  closely  approximates  akl35  is 
designated  as  the  preferred  model.  This  piece-wise  optimization  of  the  travel-time  curve  markedly  improves 
predictability,  especially  at  regional  and  teleseismic  distances. 

Figure  3  shows  a  cartoon  illustrating  our  preferred  regionalized  models  at  station  AAE.  P- wave  velocity  versus 
depth  profiles  for  each  distance  range  at  each  station  are  shown  in  Figure  4.  The  preferred  one-dimensional  regional 
models  are  then  merged  into  one  two-dimensional  model  over  the  entire  distance  rage  of  the  model,  from  regional 
out  to  teleseismic  distances  (Figure  2).  The  profiles  in  Figure  4  all  display  some  distance  and  depth  discontinuity 
artifacts  between  the  individual  models,  illustrating  the  importance  of  using  indicator  functions  when  merging  the 
regionalized  models  into  one  2-D  model. 

To  test  the  predictive  power  of  the  2-D  models,  we  compute  the  rms  and  mean  residuals  between  the  declustered  P- 
phase  arrivals  at  each  station  and  the  arrivals  predicted  by  both  the  optimized  2-D  model  and  the  akl35  model. 

Those  results  are  shown  in  Tables  2a  and  2b  below. 


Table  2a:  rms  residuals  (seconds  squared) 


STATION 

data  vs  akl35 

data  vs  best  2-D  model 

Reduction 

AAE 

2.67838 

2.05136 

23% 

ANTO 

2.16813 

1.68400 

22% 

AQU 

2.13796 

1.76301 

18% 

TBT 

2.76454 

1.69335 

39% 

Table  2b: 

mean  residuals  (seconds) 

STATION 

data  vs  akl35 

data  vs  best  2-D  model 

Improvement 

AAE 

2.12126 

0.86341 

59% 

ANTO 

1.18942 

-0.00094 

92% 

AQU 

1.22124 

-0.02453 

96% 

TBT 

0.99596 

-0.25633 

74% 

Histograms  of  the  distribution  of  residuals  are  shown  in  Figure  5.  Table  2  and  Figure  5  demonstrate  a  clear 
reduction  in  travel-time  variance  between  akl35  and  our  preferred  2-D  models.  In  addition,  our  2-D  models  provide 
a  mean  residual  closer  to  0.  These  results  suggest  that  we  are  improving  predictability  with  these  models. 


CONCLUSIONS  AND  RECOMMENDATIONS 


Summary 

An  accurate  velocity  model  of  the  earth  is  a  fundamental  component  of  seismic  event  location,  which  in  turn  is 
integral  to  the  effective  monitoring  of  the  Comprehensive  Nuclear-Test-Ban  Treaty  (CTBT).  Because  global 
velocity  models  can  be  inadequate  for  prediction  of  travel  times  at  regional  and  near  teleseismic  distances,  we  are 
developing  two-dimensional,  station-specific  velocity  models  that  are  optimized  to  predict  travel  times  for  IMS 
stations  in  the  Middle  East  and  North  Africa. 

We  develop  an  adaptive  grid  search  method  that  efficiently  samples  the  space  of  reasonable  P-wave  velocity  models 
of  the  earth.  Optimized  1-D  regionalized  models  for  each  IMS  station  are  better  able  to  predict  travel-time  arrivals 
than  akl35.  Regionalized  models  are  being  mathematically  combined  into  one  2-D  model,  using  indicator  functions 
and  smoother  methodologies  to  reduce  distance  and  depth  discontinuity  artifacts  between  individual  models.  We 
find  that  this  piece-wise  optimization  of  the  travel-time  curve  improves  predictability,  especially  at  regional  and  near 
teleseismic  distances. 

Further  research 

These  are  preliminary  results  from  an  initial  set  of  IMS  seismic  stations;  we  are  in  the  process  of  incorporating  the 
remainder  of  the  IMS  stations  into  our  study.  A  station-specific  analysis  of  the  capability  of  our  2-D  models  to 
predict  the  observations  will  help  us  identify  regions  where  the  model  needs  to  be  improved  or  alternative  models 
need  to  be  applied.  We  will  further  optimize  our  regionalized  models  by  performing  a  finer  adaptive  grid  search 
centered  near  the  residual  lows  identified  in  our  initial  analysis.  We  will  also  experiment  with  those  additional 
model  parameters  we  have  identified  as  mechanisms  to  fine  tune  the  model,  such  as  positive  and  negative  velocity 
gradients  in  the  upper  mantle. 


We  will  test  the  predictive  capability  of  each  2-D  model  by  cross-validating  the  data.  Our  cross-validation  data  set 
at  each  station  consists  of  10%  of  the  groomed  modeling  data  set.  We  expect  that  cross-validating  the  travel-time 
predictions  by  using  arrivals  from  this  set  of  independently  located  events  will  demonstrate  a  marked  reduction  in 
the  variance  of  the  travel-time  model  error  distributions. 

A  subset  of  events  from  the  1991  Racha  earthquake  sequence  recorded  at  stations  KAS,  ARU,  SVE,  KVT,  GAR ,  and 
KHO  will  be  used  to  demonstrate  improvement  in  earthquake  location  using  our  2-D  models  over  our  1-D  models 
and  over  a  generalized  earth  model  such  as  akl35. 

Finally,  we  will  use  our  travel-time  models  at  each  station  in  tandem  with  nonstationary  spatial  corrections 
(nonstationary  Bayesian  kriging)  to  further  improve  our  capability  to  accurately  locate  all  seismic  events  in  this 
region. 

This  work  was  performed  under  the  auspices  of  the  U.S.  Department  of  Energy  by  University  of  California 
Lawrence  Livermore  National  Laboratory  under  contract  No.  W-7405-Eng-48  for  the  Office  of  Research  and 
Development,  NN-20,  within  the  Office  of  Nonproliferation  and  National  Security,  NN-1.  LLNL  contribution 
UCRL-JC- 138983. 

REFERENCES 


Buland,  R.  and  C.  Chapman.  The  computation  of  seismic  travel  times.  Bull.  Seism.  Soc.  Am.,  73,  1271-1302,  1983. 

Engdahl,  E.R.,  R.  van  der  Hilst,  and  R.  Buland,  Global  teleseismic  earthquake  relocation  with  improved  travel  times 
and  procedures  for  depth  determination.  Bull.  Seism .  Soc.  Am.,  88,  722-743,  1998. 

Kennett,  B.L.N.,  E.R.  Engdahl,  and  R.  Buland.  Constraint  on  seismic  velocities  in  the  earth  from  travel-times, 
Geophys.  J.  Int.,  122,  108-124,  1995. 


Degress  Urtmacte 


0  uqn:  L:  ■  Lu.".idt  uda 


Figure  1:  Location  of  the  27  International  Monitoring  System  (IMS)  seismic  stations  in  the  Middle  East  and  North 
Africa  for  which  we  have  developed  2-D  travel-time  models.  We  will  present  detailed  modeling  results  from  the 
four  IMS  stations  highlighted  by  the  encircled  triangles. 
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Figure  2:  Cartoon  depicting  the  merging  of  regionalized  1-D  models  into  a  2-D  model.  Merging  is  accomplished 
using  indicator  functions  to  reduce  distance  and  depth  discontinuity  artifacts  between  the  individual  models.  The 
result  is  a  radially  heterogeneous  and  azimuthally  invariant  (see  Figure  3)  travel-time  model  of  the  crust  and  upper 
mantle.  This  methodology  provides  optimal  models  for  the  three  distinct  ray-bottoming  depths,  allowing  increased 
predictability  and  a  smooth  travel-time  curve. 


Regionalized  radially  heterogeneous 
and  azimythally  invariant 
travel-time  models 


Example:  Station  AAE 


Regional  distance  model 

ennui  tftisfcntB*  »  35  tm 
Vp  -truvCLil  Isyer  1  •  5,5  km/ts 
Vp  •  -crusEal  laytr  £  w  6.5  km/t 
Pn  u4lacit/  s  6.0  km/'s 


Upper  mn  ill  I  c  distance  model 
crustal  thickness  ■*  55  k  fn 
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Figure  3:  The  map  in  the  upper-right-hand  corner  shows  the  azimuthal  distribution  of  a  typical  modeling  data  set. 
This  particular  data  set  was  recorded  at  station  AAE.  The  data  have  been  parsed  into  three  earthquake-station 
distance  ranges:  regional  (white  circles;  1°-13°),  upper  mantle  (gray  circles;  13°-30°),  and  teleseismic  (30°-90°). 
This  figure  also  illustrates  the  best  regionalized  radially  heterogeneous  and  azimuthally  invariant  travel-time  models 
at  the  IMS  seismic  station  AAE. 
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Figure  4:  Preferred  regionalized  models  for  stations  AAE,  ANTO,  AQU,  and  TBT.  The  P-wave  velocity 
versus  depth  profiles  for  each  distance  range  at  each  station  are  shown.  The  solid  line  corresponds  to  the 
regional -distance  model;  finely  dashed  line  corresponds  to  the  upper  mantle-distance  model;  coarsely 
dashed  line  corresponds  to  the  teleseismic-distance  model. 


bumtoet  or  a&tinanco<i  Hu  isto  ■  or  iftxrJiranooe  ■- jurer-  or  aoanwioee  Hlmflef  or  oeaKreflcee 


Station  AAE 


Station  ANTO 


ll^ll^■lllll^l .  ...  :■:■•  ■iimr 

It  I  t  «  1  1  1  4  £  B  !■- 

TtivsHImt  vrftdtjiil 


1(2  :  b  ■!  I  J  a  ■■:  «  1  1L- 


7  raiwHhwe  raalcfciaf  (Mconcts). 


TraitsHUde  reaidjat  i  seconds! 


TtimHIma  r..-»Judl  (seoafida) 


Station  AQU 


Fraust-tima  r-saitUai  iaaodfitfet 


TrsveMice  nshidufr  •a&odiv.tet 


Stsion  TBT 


:>i  £  4  ^  3  5  I  4  4  B  1 U 


FraveHifne  r-iisitua^Krdit^st 


Figure  5:  Distribution  of  travel-time  residuals  predicted  by  a  global  velocity  model  ( akl35 )  and  2-D 
models  computed  for  IMS  seismic  stations  AAE ,  ANTO ,  AQU,  and  TBT.  Rrms  =  rms  residual  (seconds 
squared);  Rmean  =  mean  residual  (seconds). 


